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ADVANCES  IN  THE  NUMERICAL  COMPUTATION 
OF  CAPILLARY-GRAVITY  WAVES 


Jean-Marc  Vanden-Broeck 

Stanford  University 
Stanford,  California 


§1.  INTRODUCTION 

This  paper  deals  with  the  computation  of  symmetric  finite 
amplitude  waves  propagating  without  change  of  form  on  the  sur¬ 
face  of  a  liquid  above  a  horizontal  flat  bottom.  VJe  assume 
the  liquid  to  be  inviscid  and  incompressible  and  the  flow  to 
be  irrotational .  The  free-surface  condition,  including  the 
effects  of  capillarity,  is  used  in  its  exact  nonlinear  form. 

In  Sec.  2  we  formulate  the  problem  as  an  integro-dif f er- 
ential  equation  system  for  the  unknown  shape  of  the  free  sur¬ 
face.  This  system  consists  of  a  nonlinear  differential  equa¬ 
tion  coupled  with  a  linear  integral  equation.  A  numerical 
scheme  based  on  Newton's  iterations  is  derived  to  solve  these 
equations.  Details  of  the  numerical  procedure  are  given  in 
Sec.  3.  The  formulation  of  the  problem  and  the  numerical 
method  used  to  solve  it,  follows  closely  the  work  of  Schwartz 
and  Vanden-Broeck  [8]  and  Vanden-Broeck  and  Schwartz  [10]. 

In  recent  years  important  progress  has  been  achieved  in 
the  calculation  of  steep  gravity  waves  in  water  of  arbitrary 
uniform  depth.  For  example,  Schwartz  [7]  extended  Stokes' 
series  to  high  order  by  computer  methods  and  then  recast  these 
polynomials  as  Pad6  approximants .  High  accuracy  solutions 
were  obtained  in  that  way.  This  a;  ^roach  was  further  refined 
by  Longuet-Higgins  [4]  in  the  infinite  depth  case  and  by 
Cokelet  [2]  in  the  finite  depth  case.  In  Sec.  4.1  we  use  the 
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numerical  procedure  of  Secs.  2  and  3  to  compute  steep  gravity 
waves  in  shallow  water.  These  computations  numerically  con¬ 
firm  the  validity  of  the  use  of  Pade  approximants  as  applied 
to  gravity  waves.  While  the  convergence  of  Coke let's  series 
deteriorates  for  steep  waves  in  very  shallow  water,  our  numer¬ 
ical  scheme  remains  efficient  for  depths  as  small  as  1/120  of 
a  wavelength. 

Approximate  solutions  for  gravity-capillary  waves  were 
published  long  ago  by  Harrison  [3].  The  surface  profile  was 
sought  as  a  Fourier  series  in  the  horizontal  coordinate  with 
coefficients  that  are  power  series  in  the  wave  amplitude. 

This  perturbation  expansion  invoked  Stokes'  hypothesis  that 
the  n-th  Fourier  coefficient  is  of  n-th  order  in  the  amplitude. 
The  use  of  this  hypothesis  resulted  in  infinite  values  for  cer¬ 
tain  of  the  series  coefficients  when  the  dimensionless  capil¬ 
lary  number  k  is  the  reciprocal  of  an  integer  greater  than  one. 
For  the  particular  case  k  =  1/2,  Wilton  [11],  by  revoking 
Stokes'  hypothesis  and  reordering  the  terms  of  the  series,  was 
able  to  find  two  solutions.  This  analytical  approach  was  fur¬ 
ther  extended  by  Pierson  and  Fife  [6],  Nayfeh  [5],  and  Chen 
and  Saffman  [1].  The  numerical  scheme  of  Sec.  3  was  applied 
by  Schwartz  and  Vanden-Broeck  [8]  to  compute  gravity-capillary 
waves  in  deep  water.  A  summary  of  their  results  is  given  in 
Sec.  4.2.  In  addition,  we  show  that  the  wave  speed  parameter 
p  is  not  a  single  valued  function  of  k.  This  fact  enables  us 
to  extend  Schwartz  and  Vanden-Broeck ' s  results. 

All  the  capillary-gravity  waves  obtained  by  Schwartz  and 
Vanden-Broeck  [8]  are  ultimately  limited  by  contact  with  adja¬ 
cent  waves.  Vanden-Broeck  and  Keller  [9]  have  developed  a 
numerical  procedure  to  construct  waves  of  higher  amplitude. 

A  discussion  of  their  method  is  given  in  Sec.  4.3. 

52.  MATHEMATICAL  FORMULATION 

We  consider  two-dimensional,  periodic  waves  of  wavelength 
X  and  phase  velocity  c  propagating  on  the  surface  of  a  liquid 
under  the  combined  effects  of  gravity  g  and  surface  tension  T 
over  a  horizontal  bottom.  We  choose  a  frame  of  reference  in 
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which  the  waves  are  steady,  as  is  the  fluid  notion,  which  is 

assumed  to  be  a  potential  flow. 

The  variables  are  made  dimensionless  by  referring  them  to 

1/2 

the  velocity  scale  (g\/2v)  7  and  to  the  length  A/2tt.  In  addi¬ 
tion,  we  introduce  the  wave-speed  parameter 


(1) 


and  the  dimensionless  capillary  number 


K 


4tt2T 

pgA 


(2) 


Here  p  is  the  fluid  density. 

The  condition  of  constant  pressure  (p  =  0)  on  the  free 
surface  can  be  written 


i  qq*  +  V  +  |  =  f  (3) 

Here,  q  =  u  -  iv  is  the  complex  velocity  and  the  asterisk  sig¬ 
nifies  complex  conjugation.  The  acceleration  of  gravity  acts 
in  the  negative  y  direction.  The  y-axis  is  chosen  as  a  line 
of  symmetry  of  a  surface  wave.  R  is  the  radius  of  curvature 
counted  positive  when  the  center  of  curvature  lies  on  the 
fluid  side  of  the  free  surface. 

The  choice  of  the  Bernoulli  constant  in  (3)  fixes  the 
origin  of  y  as  the  undisturbed  level  of  the  free  surface  for 
which  the  velocity  and  the  curvature  are,  respectively,  equal 
to  and  zero. 

Let  the  stream  function  assume  the  values  zero  and  -Q  on 
the  free  surface  and  on  the  bottom,  respectively.  The  undis¬ 
turbed  fluid  depth  d  is  defined  by 

d  =  -5-  (4) 

We  choose  the  complex  potential 
f  =  4>  +  i'/' 

as  the  independent  variable. 

In  order  to  satisfy  the  boundary  condition  3y / 3 tp  =  0  on 
the  bottom  =  -Q,  we  reflect  the  flow  in  the  boundary  ip  =  -Q. 
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Thus  we  seek  z  =  x  +  ty  as  an  analytic  function  of  f  in  the 
strip  -2Q  5  ip  s  0. 

It  is  convenient  to  introduce  the  following  change  of 
variables 

f  =  <j>  +  iiJ;=i/ij  log  £  (5) 

where  z;  =  rei9.  Relation  (5)  maps  the  bottom  i/;  =  -Q,  the  free 

surface  ip  =  0,  and  its  image  ^  =  -2Q,  respectively,  onto  the 

circles  r  =  rg  =  e  r=l,  and  r  =  r^.  The  values  of  the 

real  and  imaginary  parts  of  the  function  F(t)  =  tdz/dc  -  i  on 

2 

the  free  surface  r  =  1  and  on  its  image  r  =  r^  are  related  by 
the  ■•'.entities 


x'(0)  =  -1  -  IF  (el8 )  =  -1  -  IF  (r^e10 )  (6) 

y ' (9 )  =  RF(ei9)  =  -RF  (r^e19 )  (7) 

Here,  x'(0)  and  y'(9)  represent,  respectively,  the  real  and 
imaginary  parts  of  the  function  dz/d0  on  the  free  surface  r  =  1. 
These  variables  have  the  periodicity  property 

r+V  r+TT 

I  [X'  (6)  +  1]  de  -  y'  (6)  d6  =  0  (8) 

'  -TT  I  -IT 


In  order  to  find  a  relation  between  x’ (6)  and  y’(0)  we 

apply  Cauchy's  theorem  to  the  function  F(c)  in  the  annulus 
2 

Tq  5  |?|  s  1.  Using  the  relations  (6),  (7),  and  (8)  we  find, 
after  some  algebra, 


x'  (9)  +  i  =  27  [  y’  «.)  cot  ~-9)  d4> 

*  -TT 

r q  r tt  [1  +  x'  (0)]  [rjj  -  cos(0  -  <J> )  ]  - 
11  J  -Tt  r^  -  2r^  cos  (0  -  <}>)” 


y '  ( <J> )  sin  (<t>  -  6) 

+  1 


d<t> 


(9) 


the  first  integral  being  of  Cauchy  principal  value  form. 
Exploiting  the  bilateral  symmetry  of  the  wave  about  0  =  0,  we 
rewrite  (9)  as 
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if11  1  1 

x'(6)  +  1  =  jjj-  J  ^  y '  ($ >  (cot  |(4  -  0)  +  cot  +  6)  )d$ 

Tg  ,n  (1  +  x'(4)][r2  -  cos(0  -  4)]  -  y '  (0)  sin  (<(i  -  6) 

n  '  0  r0  ”  2r0  cos^6  “  <l>)  +  1 

r2  ,t\  (1  +  x'  (4)  ]  [r2  -  cos(6  +  0)]  -  y'(<t>)sin(4 
—  — —  - - - - - - - 

•*0  r.  -  2r„  cos  (0  +  41 )  +  1 


J0  r0  -  2rQ  cos  (6  -  4)  +  1 

r2  /-tt  [1  +  x'  (4)  ]  [r2  -  cos(6  +  4)]  -  y'(4)sin(4  4  6) 


rQ  -  2rQ  cos (6  +  4)  +  1 


The  surface  condition  (3)  can  now  be  rewritten  as 
y  +  §  f(x'2  +  y’2)-1  -  1] 

+  k (x'y"  -  y'x") (x1 2  +  y,2)"3/2  =  0 


In  addition  to  the  parameters  r^,  k,  and  y,  a  given  wave 
is  characterized  by  a  third  parameter  s,  which  is  a  measure  of 
the  wave  steepness.  This  parameter  can  be  defined  in  many 
different  ways.  Precise  definitions  appropriate  to  each  of 
the  problems  considered  will  be  given  in  the  following  sec¬ 
tions.  Dimensional  analysis  implies  that  a  functional  rela¬ 
tionship  should  exist  among  these  parameters: 

f(y,K,s,rQ)  =  0 

We  shall  consider  two  closely  related  numerical  schemes. 
In  the  numerical  scheme  I,  we  fix  r^,  s,  and  <  and  we  seek  the 
function  y(0),  0  S  [0,ir]  and  a  value  for  y.  In  the  numerical 
scheme  II,  we  fix  r^,  s,  and  y  and  we  seek  the  function  y(6) 
and  a  value  for  k. 


§3.  NUMERICAL  PROCEDURE 

We  seek  a  numerical  solution  of  the  integro-dif f erential 
system  of  Eqs.  (10)  and  (11)  by  a  finite  difference  method. 
Introducing  a  uniform  mesh,  we  have 

6i  =  ((i  -  1)/N]tt  i  =  1,  ...,  N  +  1  (12) 


Since  the  wave  is  symmetrical,  we  have 
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Thus  the  unknown  function  3y/30  can  be  represented  by  the  vec¬ 
tor  of  dimension  N  -  1, 

y'  =  (yfi  »  •  ••*  ya  )  (13) 

®2  0N 


where 


ye 


0=0, 


We  now  define  the  midpoints  y,  =  i  (0,  +  9,+,),  i  =  1, 

N,  and  we  represent  the  values  of  3x/39,  3y/30,  and  y  at 
the  points  y ,  by  the  vectors  x^,  y^,  and  y  . 

We  seek  to  satisfy  the  system  of  Eqs.  (10)  and  (11)  at 
the  N  points  y,.  The  integrals  in  (10)  are  evaluated  at  the 
points  y,  by  the  trapezoidal  rule  (which  is  of  infinite  order 
since  the  integrand  is  periodic) .  The  integrals  involving 
the  function  y'(<f>)  are  computed  by  integration  over  0,.  The 
singularity  of  the  Cauchy  principal  value  is  automatically 
taken  into  account  since  the  quadrature  is  symmetrical  with 
respect  to  the  singularity.  The  remaining  integrals  in  (10) 
are  computed  by  integration  over  y, .  Thus,  we  obtain 

1  +  x^  =  Ay'  +  B(x^  +  1) 


or 

XjJj  =  —1  +  (I  -  B)  _1Ay '  (14) 

Here,  A  and  B  are  known  matrices  and  I  is  the  unit  matrix.  We 
note  that  B  =  0  for  infinite  depth  (i.e.,  r^  =  0)  so  that  no 
matrix  inversion  is  needed  in  this  particular  case. 

The  vectors  y^  and  yM  are  expressed  in  terms  of  y '  by  a 
sixth-order  interpolation  formula  and  by  a  sixth-order  quadra¬ 
ture  formula,  thus, 


In  -  cr 

(15) 

Xm  =  yo  +  Dy' 

(16) 

where  C  and  D  are  known  matrices.  The  elevation  y^  of  the 
free  surface  at  0  =  0  has  to  be  found  as  part  of  the  solution. 
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Substituting  (14),  (15),  and  (16)  into  (11),  we  obtain  a  sys¬ 

tem  of  N  nonlinear  algebraic  equations.  Using  the  definition 
of  s  we  obtain  an  extra  equation.  Thus  we  have  N  +  1  equa¬ 
tions  for  the  N  +  1  unknowns  (y^,  y^,  y^,  p)  in  scheme  I 

and  (y^ ,  ...,  y^,  yQ,  k)  in  scheme  II.  This  system  of  equa¬ 

tions  was  solved  by  Newton's  iterations.  Each  iteration  re¬ 
quires  the  inversion  of  a  matrix.  We  note  that  the  matrices 

A  and  B  in  (14)  depend  only  on  rn.  Thus  computing  time  can  be 

_1  u 

saved  by  computing  (I  -  B)  A  at  the  beginning  of  the  program. 

In  all  the  cases  to  be  discussed,  the  numerical  scheme 

was  found  to  converge  quadratically .  Four  or  five  iterations 

were  required  to  satisfy  the  alaebraic  equations  with  an  error 
-12 

less  than  10  .  With  N  =  80,  each  iteration  takes  about  3.5 

seconds  on  a  CDC  6600  computer. 


§4.  DISCUSSION  OF  RESULTS 


Gravity  Waves  in  Shallow  Water 


Vanden-Broeck  and  Schwartz  [10]  have  applied  the  numeri¬ 
cal  scheme  of  Sec.  3  to  compute  steep  gravity  waves  in  shallow 
water.  The  effect  of  surface  tension  was  neglected  li.e., 
k  =  0)  and  the  parameter  s  was  chosen  to  be 


(17) 


Here,  yc  is  the  elevation  of  the  crests  of  the  wave.  For  the 
highest  wave,  the  velocity  at  the  crest  vanishes.  Thus,  from 
(3),  yc  =  u/ 2,  so  s  =  1  for  the  highest  wave.  In  general  s 
ranges  between  0  and  1. 

In  order  to  compute  accurately  the  crests  of  the  waves, 
which  become  sharp  when  s  increases,  a  new  independent  varia¬ 
ble  8  was  introduced  by  the  relation 

0  =  8  -  a  sin  8  (18) 

The  closer  a  is  to  one,  the  greater  the  concentration  of  mesh 
points  near  the  crest.  It  was  found  that  steep  waves  could  be 
computed  by  choosing  a  =  0.999. 

Table  1  shows  values  of  u  for  r^  =  0.5  computed  with  40, 
60,  and  80  mesh  points.  The  computed  values  for  N  =  80  have 
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TABLE  1. 

Values  of  p 

for  r^  =  0.5 

and  0.53  s  s  5 

0.992. 

s 

N  =  40 

N  =  60 

N  =  80 

Cokelet 

0.52966 

0.666501 

0.666501 

0.666501 

0.666501 

0.69331 

0.706443 

0.706443 

0.706443 

0.706443 

0.84727 

0.748230 

0.748230 

0.748230 

0.748230 

0.92326 

0.764404 

0.764403 

0.764403 

0.764403 

0.96149 

0.767749 

0.767750 

0.767750 

0.767748 

0.97687 

0.767537 

0.767540 

0.767540 

0.76754 

0.98458 

0.767089 

0.767096 

0.767097 

0.76707 

0.99228 

0.766551 

0.766556 

0.766557 

0.76648 

converged  to  six  places  following  the  decimal  point  for  s  <0.97 
and  to  five  places  for  s  >  0.97.  The  last  column  contains  the 
values  of  p  obtained  by  Cokelet  [2] .  The  lowest  number  corre¬ 
sponds  to  the  highest  steepness  for  which  Cokelet  computed  the 
wave  speed.  His  results  are  in  good  agreement  with  our  numer¬ 
ical  values.  Thus,  the  validity  of  the  Pade  approximant  method 
for  finite  depth  gravity  waves  is  confirmed.  It  should  be 
pointed  out  that  Cokelet  used  another  parameter  instead  of  our 
parameter  s.  Thus,  a  part  of  the  small  discrepancy  at  the  bot¬ 
tom  of  the  table  may  be  attributed  to  the  loss  of  accuracy  in¬ 
volved  in  the  calculation  of  s  from  Cokelet 's  results. 

The  convergence  of  Cokelet' s  results  deteriorates  for 
steep  waves  in  shallow  depth.  On  the  other  hand,  our  numerical 
scheme  remains  efficient  for  depths  as  small  as  1/120  of  a 


TABLE  2. 

Values  of  p 

for  rQ  =  0.9. 

S 

N  =  60 

N  =  80 

Cokelet 

0.42183 

0.127329 

0.127329 

0.12733 

0.61793 

0.142055 

0.142054 

0.141983 

0.81142 

0.158013 

0.158011 

0.1578 

0.90849 

0.154092 

0.164091 

0.1638 

0.95 

0.165039 

0.165038 

0.164 

0.98 

0.164678 

0.164684 

0.163 

0.99 

0.164433 

0.164437 

0.162 
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wavelength.  In  Table  2  we  compare  the  two  methods  for  the 
smallest  depth  r^  =  0.9  considered  by  Cokelet.  For  steep 
waves  our  results  have  converged  to  five  decimal  places. 
Cokelet' s  series  expansion  method  gives  only  two  correct  deci¬ 
mal  places  for  steep  waves. 

§4.2  Gravity-Capillary  Waves  in  Deep  Water 

The  procedure  of  Sec.  3  was  used  by  Schwartz  and  Vanden- 
Broeck  [8]  to  compute  capillary-gravity  waves  in  deep  water. 

The  parameter  s  was  defined  as  the  steepness  of  the  wave  (i.e., 
the  difference  of  ordinates  between  a  crest  and  a  trough  div¬ 
ided  by  the  wavelength) .  Four  families,  each  labeled  with  a 
"type  number,"  were  studied.  The  numerical  results  indicate 
that  Harrison's  [3]  series  solution  agrees,  at  least  for  small 
finite  values  of  the  steepness,  with  the  family  of  type  n  when 
k  lies  in  the  open  interval  (1/n  +  1,  1/n)  for  n  ;  2  and  in 
the  interval  (1/2,  °°)  for  n  =  1.  Typical  profiles  of  the  free 
surface  for  various  values  of  y ,  k,  and  s  can  be  found  in 
Schwartz  and  Vanden-Broeck  (8). 

Figure  1  shows  the  variation  of  the  parameter  y  with  k 
for  s  =  0.03.  The  four  continuous  families  are  displayed. 

Also  shown  is  the  well-known  infinitesimal  wave  solution 

y  =  1  +  K 

Schwartz  and  Vanden-Broeck  [8]  started  the  iterations 
with  a  simple  cosine  profile.  The  numerical  scheme  was  found 
to  converge  to  the  family  of  type  n  when  <  f  (1/n  +  1,  1/n) 
for  n  ;  2  and  <  (  (1/2,  »)  for  n  =  1.  Once  a  given  wave  of 

type  n  had  been  obtained,  the  family  n  was  then  computed  by  a 
type  of  "boot-strap"  technique,  using  the  numerical  scheme  I. 
That  is,  a  converged  solution  for  one  value  of  <  was  used  as 
an  initial  guess  for  a  wave  with  <  altered  by  a  few  percent, 
and  so  on.  For  s  =  0.03  Schwartz  and  Vanden-Broeck  [8]  found 
that  family  3  could  not  be  computed  for  <  ■-  0.2515.  Similar¬ 
ly,  family  4  could  only  be  computed  for  <  =■  0.21. 

In  the  present  paper  we  have  extended  the  computations 
of  families  3  and  4  by  a  "boot-strap"  technique  using  the 
numerical  scheme  II.  For  example  .ne  solution  of  family  3 
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FIG.  1.  Variation  of  speed  parameter  m  with  capillary  number 
■  for  s  =  0.03.  The  solid  curves  are  results  of  Schwartz  and 
Vanden-Broeck  [8]  while  the  dashed  lines  show  the  results  of 
the  present  computations. 


for  •  =  0.2515  and  .  =  1.2478  was  used  as  an  initial  guess  to 
compute  the  solution  for  a  value  of  i.  slightly  less  and  so  on 
The  results  are  shown  in  Fig.  1.  Families  3  and  4  could  be 
extended  to  larger  ranges  than  those  presented  in  the  figure. 
However,  these  results  are  not  reported  here  since  our  pur¬ 
pose  is  simply  to  show  that  the  parameter  is  not  a  single 
valued  function  of 
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All  the  capillary-gravity  waves  computed  by  Schwartz  and 
Vanden-Broeck  [8]  are  ultimately  limited  by  contact  with  adja¬ 
cent  waves  for  some  value  s*  of  the  amplitude  parameter.  The 
analytic  continuation  of  these  solutions  for  s  >  s*  yield  over 
lapping  waves  with  multiple  valued  velocities.  Thus,  these 
solutions  are  not  admissible  as  solutions  of  the  physical  prob 
lem.  Vanden-Broeck  and  Keller  [9]  have  shown  how  to  modify 
the  numerical  scheme  of  Secs.  2  and  3  in  order  to  obtain  phys¬ 
ically  acceptable  solutions  for  s  >  s*.  The  general  idea  is 
to  prevent  the  free  surface  from  crossing  itself.  The  solu¬ 
tions  for  s  >  s*  have  adjacent  waves  touching  at  one  point, 
just  like  in  Schwartz  and  Vanden-Broeck ' s  highest  waves.  Thus 
each  pair  of  adjacent  waves  enclose  a  region  devoid  of  fluid, 
which  we  call  a  bubble.  The  pressure  in  the  bubble  is  a  func¬ 
tion  of  the  amplitude  parameter. 

Vanden-Broeck  and  Keller  [9]  have  performed  accurate  com¬ 
putations  in  the  particular  case  of  pure  capillary  waves  in 
water  of  infinite  depth.  They  have  shown  that  the  new  family 
of  solutions  exists  for  arbitrarily  large  values  of  the 
steepness . 
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